^<A NAME='T V ></A><A HREF— f #2^^S 



DOCUMENT-IDENTIFIER: W<A NAME='T i x/A><A HREF="#2~SPAN CLASl... Page 1 of 15 

US-PAT-NO: 5570321 

DOCUMENT- TTQ^7moi a 

IDENTIFIER: U & ^ m5ZV 

TITLE* Seismic velocity model optimization method using simulated annearling to 

determine prestack travel-times 



US PATENT NO.-PN(l): 
5570321 

US Document Identifier - DID (1): 

US 5570321 A 
TITLE - TI (1): 

Seismic velocity model optimization method using simulated annearling to determine prestack travel- 
times 



Abstract Text - ABTX (1): 

A method of determining a reference seismic v elo c ity model that includes lateral and vertical 
variations in a stratified subterranean medium is disclosed. The invention recognizes that each unknown 
prestack traveltime curve may be expressed as the sum of a best fit hyperbola and a non-hyperbolic 
perturbation term. After the best fit hyperbola has been determined and used to correct the traces for 
normal moveout (NMO), the non-hyperbolic perturbation term is solved for using a simulated annealing 
technique. In addition, the problem of topologically representing a subsurface having complicated 
horizon geometries is bypassed by assuming a smooth and continuous, rather than a layered and 
discontinuous, velocit y model. This assumption facilitates raytracing and in most cases, results in a 
model that more accurately represents the velocit y function. Error minimization during model updates is 
achieved by solving a linear system of equations subject to a set of constraint equations. 

Brief Summary Text - BSTX (2): 

This invention relates to an improved method of estimating seismic v e lo cities for use in processing 
and interpreting seismic data indicating valuable characteristics of subsurface earth formations, and 
particularly, to a method in which prestack traveltimes are determined by a simulated annealing 
technique and used to recover a geologically constrained, continuous (or smooth) subsurface v eloci ty 
model. 

Brief Summary Text - BSTX (4): 

Conventional seismic data acquisition techniques involve the use of an appropriate signal source to 
generate seismic energy and a set of corresponding receivers, such as geophones, spaced along or near 
the surface, to detect any reflected seismic signals due to seismic energy striking subsurface geological 
boundaries. The seismic signals are generated sequentially at each of a number of points along a seismic 
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prospecting path while reflections are recorded at all of the points following generation of each signal. 
The recorded signals are then organized into gathers of traces each corresponding to a common depth 
point or common midpoint along the prospect path. The basic purpose of this exploration technique is to 
allow the signals within each gather to be combined to improve the signal to noise ratio. Due to the 
different path lengths involved in each source and receiver pair, corrections must be made to the 
individual traces within the gather to place them in alignment before stacking. These corrections are 
performed by a processing technique referred to as hyperbolic normal moveout (NMO), the accuracy of 
which depends primarily on estimated velocities of the signals passing through the earth formations. 

Brief Summary Text - BSTX (5): 

Accurate determination of the velocit y distribution of the subsurface is necessary for obtaining 
accurate images of subsurface formations because errors in velocit y estimation result in errors in the 
alignment of these signals and thereby reduce the signal-to-noise ratio of the resulting stacked signal. 
Traditional seismic processing techniques, such as ray map migration, and new developments, such as 
amplitude versus offset (AVO) and multiparameter inversion techniques, are critically dependent upon 
the low-frequency v elo city field. Prestack reflection traveltimes of seismic signals are required input for 
any algorithm attempting to map accurately the veloci ty distribution of the subsurface. The 
nonhyperbolic nature of prestack traveltime curves even in the presence of flat bedding-plane geometry 
has long been recognized, and further, large model-dependent oscillations can be observed, especially at 
large offset-to-depth ratios. Compounding the problem, lateral veloci ty variations and data acquisition 
complications, such as streamer feathering, can introduce traveltime perturbations of arbitrary, non- 
hyperbolic shape. Although these problems are well known, for reasons of efficiency and signal-to-noise 
ratio considerations, stacking v e lo c it y analyses typically are still carried out as a two-parameter 
hyperbolic search. 

Brief Summary Text - BSTX (6): 

A recently developed strategy for determining velocit y variations in a stratified earth, commonly 
referred to as tomography, has been used to produce enhanced subsurface images. Traveltime 
tomography techniques involve three steps: identifying a number of key horizons in a stacked section; 
determining the corresponding prestacked traveltimes; and solving for a velocity -depth model that 
reproduces the observed traveltimes. 

Brief Summary Text - BSTX (7): 

However, tomography is not routinely practiced because of its requirement for accurate prestack 
reflection traveltime selection. Three different schemes are usually employed for selecting, or "picking," 
traveltimes, each of which has certain disadvantages. One approach is to digitize prestacked traveltime 
data on an interactive work station. However, this is a very tedious process requiring the display and 
selection of thousands of seismic events. A second approach is to use stacking vel o citi e s and zero-offset 
times to construct a subsurface macro model. However, this technique is based essentially on a 
hyperbolic approximation to arrival times, and is therefore inaccurate. A third approach is to overcome 
the need for accurate traveltime picking by employing a coherency inversion method which does not 
depend on prestack time picking and is not based on curve fitting or hyperbolic approximations. The 
method is formulated as one of global optimization of some energy function. An optimization algorithm 
produces a velocit y model which maximizes some measure of coherency. This measure is calculated on 
unstacked trace gathers in a time window along the traveltime curves which are determined by tracing 
rays through the model. Knowledge of zero-offset traveltimes for principal reflectors (for example, from 
post-stack picking) is used alternatively to update v elocities using the coherency measure and the 
interface position using zero offset time information until an optimal solution is obtained. Since velocity 
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updating using the coherency measure is a highly nonlinear process, it is performed using a type of 
Monte Carlo technique referred to as the generalized simulated annealing method for updating the 
velocity field. Coherency optimization by simulated annealing is described in detail in E. Landa et al, 
"Reference Velocity Model Estimation from Prestack Waveforms: Coherency Optimization by 
Simulated Annealing," 54 Geophysics 984-990 (1989). However, a major disadvantage of the foregoing 
coherency optimization method is that it requires large computational resources to perform computations 
for a large number of model parameters, with the ever present risk of inaccuracy due to cycle skipping 
and convergence to local minima. 

Brief Summary Text - BSTX (8): 

What is needed is a tomographic technique for analyzing seismic data in which prestack traveltimes 
are determined by a method other than manual picking, since prestack traveltimes are a required input 
for any algorithm attempting to map accurately the ve lo city distribution of a subsurface. Such a 
technique would enable reference velocit y model estimation based upon a traveltime-based solution as a 
preferred alternative to coherency optimization for the reasons described above. 

Brief Summary Text - BSTX (10): 

The foregoing problems are solved and a technical advance is achieved by a method of the present 
invention for optimization of a seismic ve l ocit y model. According to the method of the present 
invention, each unknown prestack traveltime curve may be expressed as the sum of a best fit hyperbola 
and a non-hyperbolic perturbation term. In a departure from the art, after the best fit hyperbola has been 
determined and used to correct the traces for normal m ov eo ut (NMO), the non-hyperbolic perturbation 
term is solved for using a simulated annealing technique. In a further departure from the art, the problem 
of topologically representing a subsurface having complicated horizon geometries is bypassed by 
assuming a smooth and continuous, rather than a layered and discontinuous, velo city model. This 
assumption facilitates raytracing and in most cases, results in a model that more accurately represents 
the velocity function. Error minimization during model updates is achieved by solving a linear system of 
equations subject to a set of constraint equations. 

Brief Summary Text - BSTX (11): 

In a preferred implementation, seismic events are selected and normal incidence traveltimes T.sub.OL 
(x) and stacking velocities V.sub.stL (x) are determined from acquired seismic data for the selected 
seismic events. For each selected event, correct ion for NMO is performed using a best- fit hyperbola. 
The residual NMO is determined using a simulated annealing technique, which is an optimization 
method in which a multi-parameter model space is sampled randomly at different points. Simulated 
annealing is performed independently on every seismic event and every CDP gather. 

Brief Summary Text - BSTX (12): 

The true prestack traveltimes T.sub.L (x, h) for each event are defined as the sum of the best-fit 
hyperbola and a perturbation term determined by simulated annealing. An initial ve l o city model is 
constructed on the assumption that the subsurface horizon being modeled is smooth and continuous. Ray 
map migration is performed to convert the T.sub.OL (x) normal incidence traveltime to depth. 
Subsequently, ray tracing is used to calculate model traveltimes Tcal(x, h). The model traveltimes are 
compared with the true prestack traveltimes and if the two are sufficiently close, the method is complete, 
because the model is sufficiently accurate. Otherwise, the velocit y model is updated according to a set of 
linear equations subject to a set of constraints so that the smooth vel oc ity model mimics the observed 
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layers on the seismic section. As a result, the method of the present invention takes into account the idea 
that prestack traveltime determination is a local, ensemble-based energy optimization problem that can 
be efficiently solved using a modified simulated annealing technique. 

Brief Summary Text - BSTX (14): 

Another technical advantage achieved with the present invention is that, because the velocit y is not 
solved for one layer at a time (layer stripping), but globally, deeper data is allowed to constrain the 
solution at shallower parts of the model, so that propagation and magnification of error commonly 
associated with layer stripping is bypassed. 

Brief Summary Text - BSTX (16): 

A further technical advantage achieved with the invention is that it significantly decreases the amount 
of t ime necessary to develop an accurate model. 

Drawing Description Text - DRTX (5): 

FIG. 3b illustrates the seismic event of FIG. 3a after correction for normal moveout (NMO) using 
prior art techniques. 

Drawing Description Text - DRTX (6): 

FIG. 3c illustrates the seismic event of FIG. 3a after co rrection for residual NMO using a simulated 
annealing technique of the present invention. 

Drawing Description Text - DRTX (7): 

FIG. 4 is a flow chart illustrating a method of the present invention for optimizing a reference seismic 
v e l o city model used to determine vel o c ity variations in a stratified subterranean structure. 

Drawing Description Text - DRTX (9): 

FIG. 6 illustrates the nine prestack seismic events of FIG. 5 from a single CDP gather after cor re cti o n 
for NMO, with each event occupying a 100 millisecond window. 

Drawing Description Text - DRTX (11): 

FIG. 8 illustrates the seismic events shown in FIG. 6 after correction for NMO using prestack 
traveltimes determined using the simulated annealing technique of the present invention. 

Drawing Description Text - DRTX (16): 

FIGS. 13a and 13b illustrate time-shifting of a non-hyperbolic seismic event based on a best fit 
hyperbolic moveout . FIG. 13b also shows the window over which the search for the nonhyperbolic 
component of the event is carried out. 

Detailed Description Text - DETX (3): 

FIG. 2 illustrates schematically the x-t seismic record produced by generating seismic energy into the 
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earth model 10 for the layer LI . Whena number of traces 12 are ordered for a common depth point 
(CDP) 14 (FIG. 1) according to the displacement (offset (h)) of the corresponding receiver to source, 
and are gathered, i.e., displayed next to each other, a characteristic dipping shape in the wavelets in the 
traces is observed. This is an artifact which occurs because the CDP traces 12 recorded at longer source- 
to-receiver offsets correspond to longer traveltimes for the seismic energy. Such a systematic shift to 
longer reflection times (t) due to increasing source-receiver offset is referred to as normal moveout or 
M NMO M . It is well known that NMO causes errors in determining both compressional and shear wave 
Yelocities which, if left uncorrected, will cause stacked amplitudes of seismic events to be misaligned, 
thereby masking the true characteristics of the reflecting interfaces (I). 

Detailed Description Text - DETX (4): 

Using conventional velocit y stacking analysis, a hyperbola 16 may be fitted to the wavelets, or signals, 
illustrated in FIG. 2, for the seismic event at the CDP 14. It is understood that the hyperbola 16 is an 
approximation only, the true seismic event being nonhyperbolic. 

Detailed Description Text - DETX (5): 

Correction for NMO is performed utilizing conventional processing techniques whereby the traces 12 
are compressed in time to compensate for the varying distance traveled by the energy from source to 
receiver. Ideally, an NMO co rre ct ion method will flatten the seismic event for each layer represented by 
the hyperbola 16, along the time (t) axis to a value representing the zero-offset traveltime (TO). 

Detailed Description Text - DETX (6): 

As illustrated in FIGS. 3a-3b, NMO correction methods in practice tend to leave some residual 
move ou t in the data, due to the non-hyperbolic nature of the seismic event. FIG. 3a shows a synthetic 
CDP record of a single seismic event. The seismic event is specified by a 25 Hertz Ricker wavelet and a 
nonhyperbolic traveltime (t) curve 300. The curve 300 represents the sum of a true hyperbola and a low 
wave number, low amplitude perturbation term. FIG. 3b shows a curve 302 representing cor rectio n of 
the seismic event for NMO using conventional techniques, leaving some residual moveout in the data. 
Specifically, a large residual moveout with a root mean square (rms) value of about 14 ms is still 
present, because of the nonhyperbolic nature of the seismic event. Uncorrected, this residual moveout 
will deteriorate the stacked trace when the traces are subsequently stacked to increase the signal-to-noise 
ratio. Furthermore, the non-hyperbolic component of the signal carries valuable seismic velocit y 
information that unless promptly recovered will remain unutilized. 

Detailed Description Text - DETX (7): 

FIG. 3c illustrates a curve 304 in which the residual NMO for the seismic event has been corrected, 
utilizing a simulated annealing technique of the present invention, discussed in detail below. The result 
is a curve 304 which is substantially flat along the true zero offset traveltime (T.sub.0) for the seismic 
event. As will be described, the invention accurately determines prestack traveltimes of selected seismic 
events by recognizing that the true prestack traveltime for each event is the sum of the best fit hyperbola 
solved by a velocity stacking analysis, and a perturbation term solved by a simulated annealing 
technique. The determination of prestack traveltimes in this fashion is easily automated. Specifically, the 
invention comprises an iterative method for optimizing a reference seismic ve l o city model in which 
prestack traveltimes of seismic events are determined using simulated annealing, described in detail 
below. 
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Detailed Description Text -DETX (8): 

FIG. 4 is a flow chart illustrating a method of the present invention for optimizing a reference seismic 
velocity model that best describes the velocity variations in the subsurface. 

Detailed Description Text - DETX (9): 

The method of the present invention commences at step 400 with the acquisition of seismic data 
through the utilization of well known seismic exploration techniques for generating compressional 
energy into a subsurface formation and recording reflected waves from subsurface reflecting interfaces. 
The data is conventionally recorded as seismograms depicting seismic traces, representing the amplitude 
of seismic reflections as a function of time (t) and distance (x) along a path of exploration in the x 
direction of the earth's surface. The traces are typically gathered into an x-t array, referred to as a shot 
record, and subsequently, the data are regathered in common midpoint gathers and stacked. The normal 
incidence traveltimes T.sub.OL (x) are picked from the stacked section and the associated stacking 
velocities V.sub.stL (x) are determined as a function of CDP surface location (x) for a number of layers 
L. T.sub.OL (x) can be easily obtained by digitizing seismic events (also referred to as "horizons") of 
interest on a stacked section. V.sub.stL (x) is determined using a conventional v elo c i ty analysis in a 
window centered around T.sub.OL (x). 

Detailed Description Text - DETX (10): 

The purpose of the velocit y analysis of step 400 is to refine the sampling of the stacking velocity field 
both in space and in time in order to maximize the signal-to-noise ratio of the samples. As with any 
Ydocity analysis technique, it involves an interpretive step of choosing among several coherence peaks. 
Preferably, the selected seismic events on the stacked section will be those with the best signal-to-noise 
ratio, which makes picking the best stacking vel o cit y an easy task. 

Detailed Description Text - DETX (11): 

FIG. 5 illustrates an example stacked section 500 from a marine data set in which nine seismic events 
502-518 of variable signal-to-noise ratio have been identified. The events 502-518 were digitized and 
the stacking velocities (V.sub.stL (x)) were determined using standard seismic processing software. 

Detailed Description Text - DETX (12): 

FIG. 6 illustrates a section 600 representing NMO corrected seismic events for a single CDP gather 
from the section 500 of FIG. 5. Shown are 100 ms windows 602-618 centered around the best fit 
hyperbolas for all nine events 502-518 (FIG. 5). To improve the visual perception, the windows 602-618 
are stacked on top of each other proceeding from the shallowest event in the window 602 to the deepest 
event in the window 618. If the events 502-518 were truly hyperbolic, and the hyperbolas were 
accurately recovered by conventional velocity analysis, then each event would appear as a flat line at the 
center of the corresponding window 602-618. The section 600 shows some residual m o veout, most 
prominent in windows 606 and 608. Several reasons account for the residual moveou t. which include the 
nonhyperbolic nature of the prestack traveltimes at large offset-to-depth ratios, and the unequal 
amplitude level of the events across offset. For events with decreasing amplitude with offset, the near 
traces enter the semblance calculation (discussed in detail below) with a higher weight than the far 
traces, and as a result most attention is given to aligning the near offsets rather than the far offsets. The 
event in window 606, for example, was found in some CDPs to exceed 30 ms of the residual moveout at 
the farthest offset. 
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Detailed Description Text - DETX (15): 

Referring again to FIG. 4, once the seismic data input parameters T.sub.OL (s) and V.sub.stL (x) have 
been determined in step 400, the v elocit y model optimization method of the invention proceeds to step 
402. In step 402, a residual NMO calculation is performed. Specifically, the true prestack traveltime 
curve (T.sub.L (x, h)) for seismic event L at a fixed CDP location x as a function of offset distance h is 
determined by summing the best fit hyperbola from the conventional velocity analysis with a 
perturbation term, i.e.,: ##EQU1## where .DELTA.T.sub.L (x, h) is the perturbation term. The 
hyperbolic term in the above equation is known based on the input parameters of 2-way normal 
incidence traveltime T.sub.OL (x) and stacking velocit y V.sub.stL (x). The perturbation term, i.e., 
.DELTA.T.sub.L (x, h), constitutes the nonhyperbolic component to be recovered from the data using 
the simulated annealing method of the present invention. 

Detailed Description Text - DETX (16): 

The curve .DELTA.T.sub.L (x, h) as a function of offset h can take any arbitrary shape, depending on 
the velocity distribution of the overlying sediments. To accommodate this condition by 
expanding .DELTA.T.sub.L (x, h) in the interval (h.sub.min,h.sub.max), where h.sub.min is defined as 
the minimum offset distance between the seismic source and any of the geophones and h.sub.max is 
defined as the maximum such distance, in terms of cubic B-spline basis functions F.sub.i (h): 
##EQU2## where M+3 corresponds to the total number of basis functions used and c.sub.iL (x) are the 
spline coefficients. Recovery of the nonhyperbolic component in the data .DELTA.T.sub.L (x, h) 
involves estimating the values of the spline coefficients C.sub.iL (x), for i=l to M+3. The number of 
spline coefficients may be adjusted as a function of the CDP location x and v e rt ical traveltime T.sub.O to 
reflect changes in the signal-to-noise ratio of the data. 

Detailed Description Text - DETX (17): 

In order to automate the procedure of determining .DELTA.TL(x, h), an energy function that can be 
measured along any proposed curve .DELTA.T.sub.L (x, h) must be defined so that a computer can 
perturb the curve .DELTA.T.sub.L (x, h) until the energy function is optimized (i.e., maximized or 
minimized). A proper choice for such energy function is critical for the success of the method. 
Accordingly, rather than using energy functions of semblance and power stack, which functions are 
commonly used in geophysics, a time-weighted power stack of the stacked trace with the prestacked 
trace is used. The normalized crosscorrelation approximately changes the waveform of the prestack 
seismic event the traveltime of which is being measured to zero phase and to unit amplitude at all offsets 
(FIG. 7). The tim e -weighted power stack of a zero phase normalized wavelet is an energy function the 
relief of which is accentuated and therefore its maximum or minimum is better defined and easier to 
detect. FIGS. 10a- 10c are contour plots of semblance, power stack and weighted power stack, 
respectively, as a function of T.sub.O and V, which are the two parameters that define the hyperbolic 
event illustrated in FIG. 9. 

Detailed Description Text - DETX (18): 

To maximize the time-weighted power stack of the data (FIG. 10c) along an arbitrary curve defined by 
M+3 spline coefficients would require a search for the maximum in (M+3)-dimensional space. 
Therefore, rather than carrying out a direct search for the maximum by evaluating the weighted power 
stack on a grid of points in this multidimensional space, an optimization method herein referred to as 
"simulated annealing" is used. Simulated annealing allows for a computationally efficient search, while 
at the same tim e achieving convergence of the solution to the desired global maximum rather than to a 
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less desirable local maximum. 

Detailed Description Text - DETX (25): 

Unlike other cooling schedules, in the proposed algorithm, the temperature does not decrease 
monotonically. Instead, any time a new value of E.sub.min or E.sub.max is found, the temperature is 
temporarily increased. This increase can be regarded as a scale adjustment, because for a fixed k, the 
ease with which the solution exits the valley of E.sub.min remains relatively constant, independent of 
the particular depth of the point E.sub.min. An additional advantage of the new criterion is that no 
evaluation of the multiparameter cost function E at different points prior to the beginning of the cooling 
schedule is necessary for the average properties to be measured and an initial value of the temperature to 
be determined. 

Detailed Description Text - DETX (26): 

As an example, consider the one-dimensional case presented in FIG. 11. Curve A has an overall 
concave up shape with a well defined global minimum G.sub.A and a minor secondary local minimum 
L.sub.A. Curve B has a sinusoidal overall shape with a global minimum G.sub.B and a distant major 
local minimum L.sub.B. Other local minima of secondary importance are also present in both curves. 
Curves A and B are defined analytically as the sum of a few sinusoids of different amplitude and phase. 
To illustrate the importance of the topography of the cost function on the performance of the simulated 
annealing algorithm, the same experiment of locating the global minimum was repeated one thousand 
(1000) times for each curve with different starting points x.sub.o but identical cooling schedules. The 
starting points for each curve are uniformly distributed in the x direction in the domain of definition of 
the two curves. Also, it would take an average of twenty (20) steps to cover the full length in x, 
assuming movement only in one direction. As shown in FIG. 1 1, the solution based on the Metropolis et 
al. criterion converged to point G.sub.A 991 times out of 1000 runs, but to G.sub.B only 583 times, the 
remainder being trapped in L.sub.B. By comparison, repeating the same experiment using the 
acceptance criterion and simulated annealing technique of the present invention, the solution converged 
to G.sub.A 1000 t imes and to G.sub.B 996 tim e s out of 1000 runs. 

Detailed Description Text - DETX (27): 

In FIG. 12, execution of the simulated annealing algorithm begins in step 1200, in which the initial 
temperature and the number of cooling steps are set to predetermined values The number of cooling 
steps is set to about 500 and has been proved adequate by experimentation. The initial temperature value 
is any arbitrary positive constant. It is set to a fraction of the cost function E the first time E is 
calculated. In step 1201, the variable k is set equal to 1. In step 1202, the variable i is set equal to 1. In 
step 1204, a new value for the coefficient c.sub.i, is obtained by c'.sub.i =c.sub.i -hdelta.c.sub.i. In step 
1206, a new energy coefficient E(c'.sub.i) is calculated as a weighted power stack along a traveltime 
curve that accounts for the new value of coefficient c.sub.i. In step 1208, a determination is made 
whether E(c'.sub.i) is greater than E(c.sub.i). If not, execution proceeds to step 1210, in which a 
determination is made whether the acceptance criterion is satisfied. If not, execution proceeds to step 
1212, and c.sub.i remains unchanged. 

Detailed Description Text - DETX (29): 

As illustrated in FIGS. 13a and 13b, to improve performance, the solution is bounded by a window .+- 
. W(h) centered around the best fit hyperbola from the velocit y analysis step. A hyperbolic window is 
usually used; however, other window shapes could also be used. The moveout at the far offset can be 
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easily specified after visually inspecting NMO-corrected records. Because of the use of a nonflat 
window, the domain of definition [-.DELTA.c.sub.i, +.DELTA.c.sub.i ] is different for each spline 
coefficient c.sub.i and v.sub.i such that: ##EQU5## where T.sub.O (h) are the spline basis functions. 

Detailed Description Text - DETX (34): 

where .gamma, is a real number. Ignoring points T.sub.L (x.sub.me, h) a new surface S'.sub.L (x, h) is 
refit and S'.sub.L (X.sub.me, h) is substituted for T.sub.L (x.sub.me ? h). In practice, this process is 
repeated, with one point being rejected at a time . This procedure is repeated independently for every 
interface L. Typically, the number of mispicked events does not exceed three percent, but the dropped 
events T.sub.L (x.sub.me,h) have a dramatic effect on the standard deviation of the surface misfit 
Err.sub.L . 

Detailed Description Text - DETX (35): 

Referring again to FIG. 4, after T.sub.L (x, h) have been determined for each event by simulated 
annealing, the method proceeds to step 404, in which an initial seismic velocity model V(x,z) is 
constructed. 

Detailed Description Text - DETX (36): 

In step 406, model traveltimes T.sub.cal (x, h) are determined by map migrating the picked events on 
the seismic section from time to depth and then ray tracing the model in a conventional fashion. 

Detailed Description Text - DETX (37): 

The calculated traveltimes T.sub.cal (x, h) from step 406 are compared with the observed traveltimes 
T.sub.OL (x, h) from step 402. In step 408, a determination is made whether a error norm is satisfied, 
i.e., whether .vertline.T.sub.OL -T.sub.cal .vertline. is less than a given number. If so, then the model 
provides an accurate representation of the subsurface and execution ends in step 410. Otherwise, 
execution proceeds to step 412, in which the ve l o city model is updated. It should be understood that ray 
tracing is possible because in its propagation, a ray follows certain physical laws that determine its 
change in direction. In turn, these physical laws are dependent upon the distribution of the velocity in a 
subsurface. This means that at any given point, it must be known where the ray is in the model, what the 
next interface in its path is, and what the velocity of the layer on the other side of the interface is. The 
above requirements introduce a need for a topological, i.e., relative position, representation of layers in 
the model. 

Detailed Description Text - DETX (38): 

In routine geophysical practice, the seismic velocity model is represented by layers of constant vMical 
and constant or smoothly varying lateral velocit y. The layer boundaries are usually chosen to correspond 
to the top of major geologic rock formations as interpreted from well logs or to major events identified 
on the seismic section. In most cases, these layers are present only over part of the section. They 
typically truncate against faults, or are intruded by salt, or they have been partially eroded away from 
previous exposure at the earth's surface (unconformities), or they thin out laterally (pinch out), or any 
combination of the above. The degree of complexity of a subsurface model is practically unlimited. 

Detailed Description Text - DETX (39): 
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Given that the potential structural complexity of the subsurface is virtually unlimited, such 
representation is not always easy. Furthermore, tracing rays in a complicated model can become 
computationally prohibitively expensive. Therefore, according to the method of the present invention, a 
smooth, rather than a discontinuous, velocity model is constructed in step 404 and the sense of layering 
is preserved via geometric constraints used when solving for the vel ocit y updates (step 412), as 
described below. This velocity field may be described in terms of two-dimensional spline coefficients: 
##EQU10## where c.sub.ij are the spline coefficients and A.sub.i and B.sub.j are basis functions in the x 
and y directions, respectively. 

Detailed Description Text - DETX (40): 

In estimating the velocity model, an iterative scheme is used, whereby starting with an initial guess 
velocit y model (step 402), the model is incrementally changed so that with each iteration, the predicted 
traveltimes better match the observed traveltimes (step 408). Mathematically this is achieved by solving 
a system of linear equations comprising two subsystems, one that seeks to minimize the residual error 
and another that implements the constraints, i.e., 

Detailed Description Text - DETX (41): 

where .delta.c=. delta.c.sub.ij is the amount of change in the two-dimensional matrix of spline 
coefficients that describes the ve l o c ity model, T.sub.OL -T.sub.cal is the difference between the real and 
the calculated traveltimes at the current iteration that is sought to be minimized in a least squares sense, 
and G is the Frechet derivative matrix ##EQU1 1## that expresses the change in the observed traveltimes 
corresponding to a change in the model parameters. The elements of the matrix G are calculated 
simultaneously with T.sub.cal during raytracing. 

Detailed Description Text - DETX (42): 

Referring to FIG. 14 the calculation of the elements of matrix F will be explained. After map 
migration of the seismic horizons from tim e to depth is completed, a unit vector field s.sup.T is defined 
that parallels the layers and constitutes a continuous representation of the layering at the current 
iteration. At layer boundaries, the direction of s.sup.T is defined as the tangent to the boundaries. 
Between layer boundaries, the x and y components of s.sup.T are obtained by linear interpolation in the 
depth direction. s.sup.N is defined as a unit vector normal to s.sup.T. 

Detailed Description Text - DETX (43): 

As shown in FIG. 15, the subsurface model is further subdivided into two domains, .OMEGA..sub.l, 
which is comprised of small bands bracketing the layer boundaries and the discontinuities in depth, 
and .OMEGA.. sub.2, which is comprised of the rest of the model. The velo c ity is required to be updated 
in such a way that the v eloc ity field remains fairly constant in the direction perpendicular to the layering 
in .OMEGA.. sub.2, whereas it undergoes maximum changes in .OMEGA..sub.l. In essence, the two- 
dimensional velocit y update field is required to resemble a surface that has different flexural rigidity 
between interfaces than in the vicinity of the interfaces. In mathematical terms, the gradient of the 
velocity updates is required to be parallel to the layering in .OMEGA.. sub.2 and perpendicular to the 
layering in .OMEGA..sub.l, which is achieved by minimizing: ##EQU12## 

Detailed Description Text - DETX (45): 

It is well known in the geophysical literature that certain sinusoidal components of the velocity field 
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cannot be resolved from the data. For example, a sinusoidal perturbation could be added laterally to the 
velocity of a layer without noticeable effect on the observed traveltimes. Since the data are not sensitive 
to such perturbations, the solution will be unstable unless properly damped. Suppression of these 
components of the velocity field is achieved by requiring that the curvature of the wave field parallel to 
the layering is minimized, i.e., ##EQU13## 

Detailed Description Text - DETX (47): 

which yields the updates to the spline coefficients that describe the v eloc ity function. In equation (6), 
the superscript T denotes matrix transpose and the weight .epsilon. expresses the relative importance of 
the data versus the constraints. 

Detailed Description Text - DETX (48): 

Once the velocity model has been updated in the manner described, execution returns to step 406, in 
which map migration and ray tracing is repeated for the updated model. This process is repeated until a 
determination is made in step 408 that the error norm has been satisfied, indicating that a satisfactory 
model has been achieved and execution ceases. 

Claims Text - CLTX (1): 

1 . A method of determining a reference seismic velocity model that includes lateral and vertical 
variations in a stratified subterranean medium, the method comprising: 

Claims Text - CLTX (2): 

displaying a stacked section of common midpoint gathers of seismic traces, said traces representing 
the amplitude of seismic reflections for seismic events as a function of time and common midpoint 
surface location along a path of exploration of said medium; 

Claims Text - CLTX (3): 

picking normal incidence traveltimes from said display and determining the associated stacking 
v e locities as a function of said location for selected seismic events from said stacked section; 

Claims Text -CLTX (4): 

determining true prestack traveltimes for each said selected event, by summing a term representing the 
best fit hyperbola for said event with a perturbation term representing the non-hyperbolic component 
thereof, said hyperbolic term being determined from said picked normal incidence travel times and said 
associated stacking veloci t i es , and said perturbation term being determined by simulated annealing; 

Claims Text - CLTX (5): 

constructing an initial reference seismic velocity model for said medium; 

Claims Text - CLTX (6): 

determining travel times from said velocit y model to be compared to said true prestack traveltimes by 
map migrating said picked normal incidence travel times fr om time to depth and then ray tracing said 
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seismic velocit y model; 
Claims Text - CLTX (7): 

comparing said travel t im e s determined from said v elocit y model with said true prestack travel ti m es, 
such that when they are substantially the same within a predetermined error norm, said velocit y model 
provides an accurate representation of said subterranean medium strata; 

Claims Text - CLTX (8): 

updating said velocity model when said predetermined error norm is not satisfied, said velocity model 
being updated by a constrained generalized least squares solution; and 

Claims Text -CLTX (9): 

returning to said step of determining travel times from said velocity model until said error norm is 
satisfied. 

Claims Text - CLTX (17): 

4. The method of claim 2 wherein said energy function comprises a t ime -weighted power stack of said 
stacked section and said prestack traveltimes. 

Claims Text - CLTX (18): 

5. The method of claim 1 wherein said constrained generalized least squares solution used for updating 
said velocity model is defined by: 

Claims Text - CLTX (19): 

where .delta.c is the amount of change in the two-dimensional matrix of said spline coefficients that 
describes said velocity model, T.sub.OL -T.sub.cal is the difference between the picked and the 
traveltimes determined from said velocity model at the current iteration that is sought to be minimized in 
a least squares sense, and G is the Frechet derivative matrix ##EQU14## that expresses the change in 
the picked traveltimes corresponding to a change in the parameters of said velocity model. 

Claims Text - CLTX (20): 

6. A method of operating a computer to determine a reference seismic velocity model that includes 
lateral and vertica l variations in a stratified subterranean medium, the method comprising: 

Claims Text -CLTX (21): 

storing seismic data acquired for selected seismic events, said data being recorded as seismograms 
depicting seismic traces representing the amplitude of seismic reflections for said events as a function of 
time (t) and common midpoint surface location (x) along a path of exploration of said medium; 

Claims Text - CLTX (24): 

picking normal incidence traveltimes (T.sub.OL (x)) from said display and determining the associated 



http://127.0.0.1:4343/eas20020919154758728.tmp?textjb^ 9/19/02 



DOCUMENT-IDENTIFIER: U^A NAME="l"x/A><A HREF="#2"^SPAN CLAi... Page 13 of 15 

• * 

stacking velocities (V.sub.stL (x)) as a function of said location (x) for selected seismic events (L) from 
said stacked section; 

Claims Text - CLTX (25): 

determining true prestack traveltimes (T.sub.L (x, h)) for each said selected event (L) as a function of 
said location (x) and an offset distance (h), by summing a term representing the best fit hyperbola for 
said event with a perturbation term representing the non-hyperbolic component thereof, said hyperbolic 
term being determined from said picked normal incidence travel time s and said associated stacking 
velocities, and said perturbation term being determined by simulated annealing; 

Claims Text - CLTX (26): 

constructing an initial reference seismic ve loc ity model for said medium; 

Claims Text - CLTX (27): 

determining travel times (T.subxal (x, h)) from said vdocity model to be compared to said true 
prestack traveltimes (T.sub.L (x ? h) by map migrating said picked normal incidence travel times from 
time to depth and then ray tracing said seismic velocit y model; 

Claims Text - CLTX (28): 

comparing said travel ti m e s (T.sub.cal (x, h)) determined from said velocity model with said true 
prestack travel times (T.sub.L (x, h)), such that when they are substantially the same within a 
predetermined error norm, said velocity model provides an accurate representation of said subterranean 
medium strata; 

Claims Text - CLTX (29): 

updating said ve l o c i ty model when said predetermined error norm is not satisfied, said velocity model 
being updated by a constrained generalized least squares solution; and 

Claims Text - CLTX (30): 

returning to said step of determining travel t imes from said v eloci ty model until said error norm is 
satisfied. 

Claims Text - CLTX (43): 

8. The method of claim 7 wherein said acceptance criteria is determined for transitions from point i to 
point j in said velocity model space that deteriorate the solution with a probability: ##EQU15## where 
p.sub.k =f(E.sub.max -E.sub.min) and E.sub.min and E.sub.max are the minimum and maximum values, 
respectively, of said energy function encountered so far since the beginning of said simulated annealing 
cooling process. 

Claims Text - CLTX (45): 

10. A computer program stored on a computer-readable medium for execution by a computer for 
determining a reference seismic velocity model that includes lateral and vertical variations in a stratified 
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subterranean medium, the computer program comprising: 

Claims Text - CLTX (46): 

instructions for storing seismic data acquired for selected seismic events, said data being recorded as 
seismograms depicting seismic traces representing the amplitude of seismic reflections for said events as 
a function of time (t) and common midpoint surface location (x) along a path of exploration of said 
medium; 

Claims Text - CLTX (49): 

instructions for picking normal incidence traveltimes (T.sub.OL (x)) from said display and 
determining the associated stacking velocities (V.sub.stL (x)) as a function of said location (x) for 
selected seismic events (L) from said stacked section; 

Claims Text - CLTX (50): 

instructions for determining true prestack traveltimes (T.sub.L (x, h)) for each said selected event (L) 
as a function of said location (x) and an offset distance (h), by summing a term representing the best fit 
hyperbola for said event with a perturbation term representing the non-hyperbolic component thereof, 
said hyperbolic term being determined from said picked normal incidence travel t i m es and said 
associated stacking ve l ocities, and said perturbation term being determined by simulated annealing; 

Claims Text -CLTX (51): 

instructions for constructing an initial reference seismic velocit y model for said medium; 

Claims Text - CLTX (52): 

instructions for determining travel ti m es (T.sub.cal (x, h)) from said veloc i t y model to be compared to 
said true prestack traveltimes (T.sub.L (x, h) by map migrating said picked normal incidence travel 
times from time to depth and then ray tracing said seismic velocit y model; 

Claims Text - CLTX (53): 

instructions for comparing said travel times (T.sub.cal (x 5 h)) determined from said velocit y model 
with said true prestack travel times (T.sub.L (x, h)), such that when they are substantially the same 
within a predetermined error norm, said vel oci ty model provides an accurate representation of said 
subterranean medium strata; 

Claims Text - CLTX (54): 

instructions for updating said velocity model when said predetermined error norm is not satisfied, said 
v eloc i t y model being updated by a constrained generalized least squares solution; and 

Claims Text - CLTX (55): 

instructions for returning to said instructions for determining travel times from said velocit y model 
until said error norm is satisfied. 
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Claims Text - CLTX (68): 

12. The computer program of claim 1 1 wherein said acceptance criteria is determined for transitions 
from point i to point j in said v elo ci t y model space that deteriorate the solution with a probability: 
##EQU17## where p.sub.k =f(E.sub.max -E.sub.min) and E.sub.min and E.sub.max are the minimum 
and maximum values, respectively, of said energy function encountered so far since the beginning of 
said simulated annealing cooling process. 
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